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We have studied the breakup and subsequent fluid flow in very thin films of partially wetting 
liquid on solid substrates, using molecular dynamics simulations. The liquid is made of short chain 
molecules interacting with Lennard- Jones interactions, and the solid is modeled as a clean crystal 
lattice whose atoms have thermal oscillations. Films below a critical thickness are found to exhibit a 
spontaneous spinodal-like instability leading to dry patches, as predicted theoretically and observed 
in some experiments. Liquid withdrawing from a dry patch collects in a moving rim whose fluid 
dynamics is only partially in agreement with earlier predictions. 

PACS Numbers: 68.45.Gd, 47.20.Ma, 64.60.My, 83.50.Lh 



The equilibrium configuration of a partially wetting 
liquid on a solid substrate is a drop, with finite contact 
angle given by Young's equation. If the liquid is initially 
placed on the substrate in the form of a uniform coating 
film, however, it is potentially unstable and may "dewet" 
the solid so as to contract into drops, or perhaps other 
shapes. This process [Upl is obviously very relevant to 
applications in materials containment and transport, and 
at the same time involves a number of challenging scien- 
tific issues, such as the nature of the instability, the fluid 
dynamics of the dewetting liquid, and the shape selec- 
tion problem for the final pattern. In this Letter, we 
describe the results of molecular dynamics (MD) simula- 
tions which address the first two issues. We first present 
new and independent confirmation of the spinodal dewet- 
ting scenario for very thin films pH , those of sub- micron 
thickness where gravity is negligible but Van der Waals 
interactions between the liquid and the solid are crucial. 
We observe unstable films whose behavior reflects sponta- 
neous rather than nucleated instabilities. We then con- 
sider the hydrodynamics of a liquid film undergoing a 
dewetting flow, and determine the time evolution of the 
film length and shape, and the issues of internal flow 
and contact angle variability. Here we find only partial 
agreement with earlier theoretical proposals. The pat- 
tern selection issue unfortunately involves the dynamics 
on larger scales than those readily addressed with MD, 
and is not considered here. Several previous papers have 
used molecular simulations to study dewetting H, but 
with an emphasis on different parameters and aspects of 
the instability, the behavior at shorter times than we con- 
sider, and no consideration of the hydrodynamics of the 
dewetting liquid. Furthermore, we disagree with some of 
their conclusions. 

The behavior of a thin liquid film on a solid depends on 
Van der Waals (VdW) forces || if the film is sufficiently 
thin, gravity otherwise, surface tension if its surface is not 
flat, and when motion occurs, viscosity. In the partially 
wetting case, the VdW forces favor a solid exposed to 
vapor rather than liquid, but for incompressible liquids 



any attempt to withdraw the liquid increases the liquid- 
vapor surface area. The competition between VdW and 
surface tension forces forces leads to an instability cri- 
terion that small sinusoidal perturbations of wavelength 
greater than A c = \J Aji^h^ja are unstable [nj . Here a 
is a microscopic length given by a 2 — (All — ^ls)/67T7, 
where the Aij (i,j = liquid L, solid S or vapor V) are 
the Hamaker constants, 7 is the LV surface tension, \x 
is the liquid viscosity, and ho is the unperturbed thick- 
ness. The linear stability calculation shows that there 
is a most unstable wavelength a/2A c with corresponding 
growth rate w max = 3a 4 7//ift,[j. Further analysis, us- 
ing weakly nonlinear theory JTj or lubrication modeling 
H may then be used to discuss the final patterns. A 
long-standing issue has been whether the origin of ob- 
served dewetting phenomena is simply the evolution of 
these "spontaneous" instabilities, analogous to spinodal 
decomposition, or whether nucleation at impurity sites 
or by contaminants is required. Recent experiments and 
previous simulations favor the former mechanism, and 
we wish to provide further evidence, by systematically 
varying the film thickness. The procedure will be to sim- 
ulate a clean substrate of width L with periodic boundary 
conditions, and films with various values of ho- The max- 
imum wavelength of any disturbance is on the order of 
magnitude of the substrate size, so an instability should 
occur only when X c {ho) is less than O(L). 

Our simulations are based on standard MD techniques 
applied to fluid flows [B 10 1, and the details are very sim- 
ilar to that used in our earlier studies of drop wetting 
pi[ . We consider a system with 18000 fluid molecules 
consisting of four atom FENE JL2[ chains, placed on 
a substrate in the form of an fcc-crystal with atoms 
tethered by linear springs to regular lattice sites. All 
atoms interact via two-body Lennard- Jones interactions, 



V(r) =4e[(r/<r) 



12 — c ij ( r / a )~ e ]i where i, j refer to the 
fluid and solid species present. The potential is cutoff 
at 2.5a for computational speed, but a z -3 tail is added 
to the solid-liquid interaction above the cutoff distance, 
corresponding to the interaction due to a half-space of 



solid. The temperature is 1.0, maintained by a Nose- 
Hoover thermostat, and the fluid density is 0.8a -3 . The 
solid has 34848 atoms, in the shape of a square sub- 
strate of side 102.6cr, and two fee cells (four staggered 
layers) in thickness. We prefer to use molecular chains 
rather than a monatomic liquid so as to reduce the vapor 
pressure and have a sharper interface. FENE molecules 
have non-Newtonian behavior under high-stress condi- 
tions, but here the stress is usually weak (but see below) 
and any other potential confining atoms into molecules 
would suffice. Note that the r~ 6 term in the potential 
and the z~ 3 tail are the microscopic basis of the semi- 
macroscopic VdW force, and the wettability of the solid 
is controlled by the Lennard- Jones solid-liquid interac- 
tion coefficient cls- By simulating the spreading of a 
liquid drop on this substrate, we find that cls = 1 gives 
complete wetting, while cls = 0.75 gives partial wetting. 
We observe in these microscopic simulations that the wet- 
ting and dewetting behavior is essentially the same if we 
remove the z~ 3 tail but increase cls- Similarly, although 
we prefer the use of a solid made of tethered molecules 
which can realistically exchange energy with the liquid, 
the dewetting behavior is not significantly modified if 
fixed potential sites are used. The in-plane variation of 
the wall interaction is crucial, however, because if the 
potential is purely a function of height, the liquid readily 
slips along the substrate. 
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FIG. 1. Initial stages of dewetting for various thick- 
nesses. Clockwise from upper left, the initial thicknesses 
and elapsed times are, respectively, (6.45a,300r), (5.7<t,250t), 
(5.20<t,100t), (4.95<t,25t). 

We have simulated the free evolution of uniform liq- 
uid films of various initial thicknesses, and Fig. 1 shows 
snapshots for four choices near the onset of instabil- 
ity. Each frame is a top view of a liquid film, where 
a molecule is displayed as the three-segment broken line 



joining its four atomic centers, so that darkness increases 
with thickness. The initial film is slightly thicker than 
those shown, h Q = 6 .7a, and is stable over a time 1200r, 
where r = (Ty/e/m is a characteristic molecular time (in 
the picosecond range). Thinner films, produced by skim- 
ming molecules off the surface of the liquid, are unstable; 
the different frames of the figure show films at an early 
stage of evolution just after one or more holes appear, 
corresponding to dry patches on the substrate. The solid 
is flat, aside from thermal fluctuations in its atomic po- 
sitions, and the liquid is homogeneous, so there is no nu- 
cleation. Furthermore, the locations of the dry spots for 
various films are uncorrelated. Once a dry spot appears, 
it invariably continues to grow, and eventually the liquid 
withdraws into a single round drop. Liu et al. |5j claim to 
the contrary to observe a nucleation mechanism, wherein 
dry patches smaller than some critical size heal. We find 
instead that once a dry region larger than a couple of 
molecular diameters appears, which is to say once the 
molecules on the two sides are out of mutual interaction 
range, it continues to grow. 

The effect of varying the thickness in the unstable cases 
is that decreasing ho decreases the time before a dry 
spot appears, and also decreases the spacing between 
the spots. The stability calculation summarized above 
predicts that these two quantities would vary as Hq and 
/ig, respectively. It would be prohibitive in computa- 
tion time to test these results quantitatively, but the ex- 
pected trend is clearly reproduced. A further difficulty in 
making a quantitative comparison is the unknown con- 
nection between the atomic interaction parameters and 
the Hamaker constants. There is a standard relation 
B for a monatomic liquid neglecting thermal motion, 
Aij = e<j 6 ir 2 piPjCij where pi is the density of species i, 
but not for molecules. If we use this expression anyway, 
along with measurements of 7 and p from earlier MD 
simulations, our results agree with the stability calcula- 
tions to within factors of 2-3. 

Now we turn to consideration of the fluid mechanics 
of the withdrawing liquid. Conventional theoretical and 
numerical studies of dewetting flows are plagued by sev- 
eral difficulties. In general, this is a moving boundary 
problem where disjoining pressure must be included to 
take account of VdW forces. Full Stokes equation calcu- 
lations are difficult due to the fine resolution that would 
be required near the solid to capture this term. Further- 
more, there is a moving contact line, with the usual ac- 
companying shear stress singularity |13J, and where one 
is completely ignorant of the appropriate dynamic con- 
tact angle. The singularity may be eliminated by use of 
a slip model for the boundary value of the fluids veloc- 
ity fl^| , but the correct angle and its possible variation 
with flow conditions is simply unknown. The computa- 
tional complexity may be significantly reduced by use of 
a depth-averaged or lubrication model M , but the stan- 
dard parameterization of disjoining pressure, as a single 
term proportional to the inverse square of the thickness, 
does not lead to an actual dewetting motion, but rather 



pins the edge of the film. This deficiency can be remedied 
by modifications to the potential which give a minimum 
at a small value of thickness, but in this case a residual 
film is left behind, not usually seen in experiment. 

The only quantitative theoretical attempt to consider 
the fluid mechanics of dewetting, by Brochard and col- 
laborators [H , is based on the strong simplifying assump- 
tions that the liquid collects into a growing rim as it 
dewets with negligible fluid motion elsewhere, that the 
rim is in the shape of an arc of a circle, and that partic- 
ular contact angles appear at the two edges of the rim. 
It is then possible to derive power laws for the speed and 
shape evolution of the rim. Of the assumptions made in 
this treatment, it is entirely reasonable that the liquid 
in the center of the film would be static until it is swept 
up by a withdrawing rim, but the particular assumptions 
made about the shape of the rim are less compelling. The 
theoretically unbiased MD simulations reported here pro- 
vide new information at microscopic resolution. 



state and collects in a rim which translates and grows, 
at least initially. As the film dewets, we occasionally ex- 
tend the region of substrate where c^s = 0.75 to keep it 
ahead of the advancing rim, while retaining the "wetting 
value" 1.0 elsewhere. In this way, the motion is purely 
from right to left, rather than inwards on both sides, and 
it is possible to study the behavior of a longer dewetted 
region. 
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FIG. 3. Time-dependence of the length (top) and rim 
height (bottom) of the dewetting film shown in Fig. 2. The 
units are a in the ordinates and r in the abscissa. 





FIG. 2. Evolution of a dewetting film, seen from the side, 
at times lOOr, 1250r, 2500r, 3750r and 5000r. 

As in the previous calculations, we begin with a solid 
substrate initially covered by a thin uniform liquid film, 
with the same density and temperature as above, but 
thickness 12er. In order to maximize the time and space 
interval over which the film evolves, we take the substrate 
to be a long narrow strip, 547cr x 17<r, so as to discourage 
instabilities in the spanwisc direction. The liquid film is 
first equilibrated with periodic boundary conditions and 
with cls—1, corresponding to complete wetting. A strip 
of liquid is removed "by hand" at one edge of the film 
while cls is reduced to a dewetting value 0.75 for the 
solid atoms in the right-hand quarter of the substrate 
where motion occurs. The film at this stage is shown 
from the side in Fig. 2a, and in subsequent frames the 
liquid at the edge of the film withdraws from the sub- 



We see that while the dewetted liquid indeed forms a 
rim which grows in size as it moves, only the outer edge 
is circular, while the inner boundary is more ramp-like. 
There is no second contact angle g where the rim meets 
the flat part of the film. The solid-liquid-vapor contact 
angle at the receding edge is seen to fluctuate with time, 
as usual at this scale, but its mean value is observed to 
be roughly 90° and is distinctly less than the static equi- 
librium value of 105° ± 5°. Hwang et al. || present plots 
of contact angle vs. time which suggest a smooth and 
systematic variation, but in our simulations (with larger 
scale and different parameters) substantial oscillation is 
seen. The position of the edge of the film as a function 
of time is shown in Fig. 3, and as predicted Q its speed 
is roughly constant. The glitch around time lOOOr is re- 
lated to the first appearance of an obvious rim, and the 
possible flattening at the largest times is probably due 
to gradual thickening of the immobile part of the film. 
Liu et al. [g| however find a different behavior at early 
times in other systems. The height and length of the 
rim were predicted theoretically to vary as \/t, at least 
for non-volatile liquids, but we observe a much slower 
growth of the rim. The rim is subject to avalanching 
along its inner side, and grows irregularly at a rate very 
roughly i 1 / 4 , Fig. 3, while the width of the rim region 



appears to stabilize. Additional, more microscopic, in- 
formation is provided by the simulation as well. Typical 
velocity fields are shown in the rim region and in the mid- 
dle of the film, respectively, in Fig. 4. There is slip in the 
very near vicinity of the contact line but not elsewhere. 
The shear stress and pressure do not show any statisti- 
cally significant variation, aside from a shear stress peak 
at the contact line which usually accompanies slip there 
E3. The flow is essentially confined to the rim region, as 
predicted, but the details of the velocity field are hard to 
resolve because the values are small and subject to small- 
system thermal fluctuations. While there is little bulk 
motion in the plateau region of the film, there is some 
motion in the vapor and also along the surface however; 
the snapshots in Fig. 2 show some wave-like excitations 
on the film surface, and there is a non-zero velocity at 
the top of the liquid plateau indicated in Fig. 4. 



ted fluid collects from the substrate into a rim moving 
with roughly constant velocity. The detailed fluid dy- 
namics of flow in the rim region and the shape evolu- 
tion of the front is unfortunately not yet predictable. A 
particularly fruitful area for further study might be the 
effects of molecular structure and interactions on dewet- 
ting dynamics. The advantage of having the forces and 
structures under explicit control may help in the design 
of controlled dewetting patterns. 
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FIG. 4. Velocity fields at the edge (top) and in the middle 
(bottom) of a receding film, averaged over the time interval 
3125-3150t and spatial bins of cross section 1.71a. The largest 
velocity vector shown is 0.2cr/r. 

Although in this Letter we have focused on the results 
of just two rather large scale simulation sequences, we 
emphasize that numerous smaller-sized runs with vari- 
ous alternative parameters exhibit essentially the same 
behavior. Among the variations we have explored in this 
way are different molecular sizes and bonding potentials, 
different substrate thicknesses, other values of the tem- 
perature, adding a (weak) gravitational force, and differ- 
ent methods for initialization. In particular, aside from 
the obvious expectation that decreasing the strength of 
the solid-fluid interaction enhances the rate of dewetting, 
we find that sufficiently thin dewetting films always fol- 
low the trends expected from the spinodal decomposition 
scenario and its linear stability analysis, and that dewet- 
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